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1. INTRODUCTION 


The aim of this paper is to present an introductory survey of mixed 
finite element methods* We shall deal first with the so-called mixed formu- 
lation of some problems arising in elasticity and hydrodynamics. Then we 
shall analyze the difficulties connected with the choice of appropriate finite 
element discretizations for a mixed formulation. Finally, we shall discuss 
some special techniques that are often helpful for solving the discretized 
problem. 

The notation of Ciarlet [17] is followed throughout. 


2. MIXED FORMULATIONS 

A precise and satisfactory definition of "mixed method" (or of "mixed 
formulation") does not exist. The term started in the engineering literature 
(Herrmann [33], [34]; Hellan [32]) in connection with the elasticity theory to 
denote methods based on the Hellinger-Reissner principle, in which both dis- 
placements and stresses were approximated simultaneously. Even among the 
mathematicians, such as Glowinski [29], Babuska [6], Crouzeix-Raviart [18], 
Johnson [36] , whose work can now be considered as pioneering, the term "mixed" 
was used only by Johnson in the context of plate bending problems. The term 
is now used in a much wider sense and has become rather vague. Here we will 
live with such vagueness, and we shall not try a new unsatisfactory defini- 
tion. Instead, we present a few cases: linear elliptic problems, Stokes 
equations for incompressible fluids, and linear elasticity problems. The 
first case will be dealt with in more detail because it is formally much 
simpler, while only a few essential points will be stressed for the other two 


cases 
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Example 1. Linear elliptic operators 

The use of mixed formulations for linear elliptic operators is rather 
recent and, as we shall see, is recommended only in some special cases. How- 
ever, its presentation is very simple and this makes it an ideal example. 
Consider the model problem: 

(2.1) div(A(x) grad u) ■ f in D C l£, 


( 2 . 2 ) 


(A(x) grad u) *n = g x on r Neu , 


(2.3) 


u " §0 on r Dir’ 


where: 

i) r Dir ^ r Neu “ 

r = 3D 

is a splitting of 

3D, 

the boundary of 

the 

domain 

D, ii) A00 is 

a smooth 

function on I), 

with 

A(x) _> a > 0 

for 

every 

x in D, iii) n_ 

is the 

unit outward normal 

to 

3D, and iv) f, 

gp 


gQ are given smooth functions in D and on ^eu’ r^ r respectively. 
Introducing for g = g 0 or g = 0 the manifold 

Hg (D) = {v| v € H^D); v = g on r Dir > , 


we can write the variational formulation of (2.1) - (2.3) as follows: 


find u € H 1 

g 


such that 


f A(x) grad u» grad v dx = -J fv dx + / 

D D “ r 


Neu 


gjvdr 


v V e hJ(d). 
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In order to derive the mixed formulation of (2.1) - (2.3) we introduce 
the variable: 


(2.4) 


£ = A(x_) grad u in D, 


so that (2.1) and (2.2) become, respectively: 


(2.5) 


div _p = f in D, 


( 2 . 6 ) 


£ ' ” ' *1 on r Neu- 


The formulation (2.3) - (2.6) is often called the mixed formulation of (2.1) - 
(2.3). Two reasonable variational formulations for (2.3) - (2.6) are now 
possible. The first one is: 

find u € H 1 (D) and p € (L 2 (D)) d such that: 

cr *■— 


(2.7) 


( 2 . 8 ) 


/ (A(x)) £ d x - / q» grad u dx = 0 V £ € (l/(D)) a , 

D D 


2/ T x\\d 


“/ P* grad v dx = / f v dx - / g,v dr V v € H?:(D)< 

D D r K7 

Neu 


In order to introduce the second variational formulation of (2.3) - (2.6) 
we define, for g - or g “ 0, the manifold: 


H (div;D) = {q| q € (L 2 (D)) d ; div q £ L 2 (D); q»n = g on T^) . 

o 
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The second variational formulation of (2.3) - (2.6) is now: 


find 


u e L 2 (D) 


and 


£ € H (div;D) 
g l 


such that: 


(2.9) / ( A(x) ) d x + / u div q d x « / g^q-ri dr v £ € H n (div;D), 

D D r n . U 

Dir 

(2.10) J v div P_ d£ =* J fv dx^ ¥ v € L^(D). 

D D 

The difference between (2.7) - (2.8) and (2.9) - (2.10) is clearly a simple 
integration by parts (or, if you prefer, a Greenes formula). However, it must 
be pointed out that the regularity required a priori for u and j3_ is some- 
how interchanged. This implies that, in discretizing (2.7) - (2.8), one has 
to use a continuous "u" and can use a discontinuous 1 V* , while in dis- 
cretizing (2.9) - (2.10) one can use a discontinuous "u" but must use a 
"j>" with divergence in L^(D) (and hence j)*ti continuous at the 
interfaces). Note also the inversion in the treatment of the boundary condi- 
tions. 

It is questionable whether (2.7) - (2.8) should be called a mixed formu- 
lation for problem (2.1) - (2.3). On the other hand, it seems acceptable to 
call the formulation (2.9) - (2.10) "mixed. " In general, the original formu- 
lation (2.1) - (2.3) is to be preferred. It is simpler, it uses just one 
variable, and many robust methods are based on this approximation. However, 
in some applications, the "auxiliary" unknown p_ defined in (2.4) is actually 
the more relevant physical variable and/or is the only information that has to 
be transferred into other equations that are coupled with (2.1) - (2.3). In 
such cases, the use of a mixed formulation might be preferred, as far as it 
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provides (as it often does) a better accuracy for _p_. In general, the formu- 
lation (2.9) - (2*10) is then used, since it deals with a smoother vector 
field _p_. It is often said that the crucial feature in the mixed approach is 
that it averages (A(x_))~* instead of A(x) • This is surely a better thing 
do, at least in one dimension, being connected with the homogenization theory; 
see for instance Babuska-Osborn [7]. However, dramatic improvements have been 
obtained by using (2*9) - (2.10) with a constant A(_x); see for instance 
Marini-Savini [40]. The true reason (if any) for the better behavior of the 
mixed formulations over the classical ones is still not completely under- 
stood. Practical experiences suggest the use of a mixed formulation for M bad 
behaved" problems, in which the variable jpOO is expected to be "smoother" 
than the variable u(x), but clearly this is not the whole story. 

Example 2. Incompressible fluids 

The Stokes equations for incompressible fluids are of the type: 

(2.11) -Au + grad p = f_ in D C if* 

(2.12) div u = 0 in D. 

Various kinds of boundary conditions can be used in connection with (2.11), 

(2.12) . For the sake of simplicity, we shall consider only the (physically 
uninteresting) Dirichlet boundary conditions 


= 3D. 


(2.13) 


u = 0 


on r 
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The natural variational formulation of (2.11) - (2.13) is: 
find u € (Hq(D))^ and p € L 2 (D) such that: 

(2.14) J grad _u:grad v dx - / p div v dx = / f»v dx ¥ v e (H^(D))^, 

(2.15) / q div u dx = 0 ¥ q € L 2 (D). 

D 

The formulation (2.14) - (2.15) had been used for years before the term "mixed 
method" came into use; however, it is recognized now that (2.14) - (2.15) 
behaves like a mixed formulation as far as the difficulties in finding good 
approximations are concerned. We shall also see that (2.14) - (2.15) easily 
falls into the same abstract framework that is commonly used for mixed 
methods. Hence, we are somehow allowed to consider (2.14) - (2.15) as a mixed 
formulation. 


Example 3. Linear elasticity problems 

For a vector valued function v(x) we define e(v) by: 


(2.16) 


_ 1 r 8v i ^ 9v j 
e ij 2 (axj 3x i 


) (i»j ■ !»••• >d), 


The linear elasticity equations are now: 


(2.17) a = E:e(u) 


d d 

(i.e., cr = £ £ E.. e (u)), 

ij r =l s=l ljrs rs “ 


div a = f 


(2.18) 


in D 
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Substituting (2.16), (2.17) into (2.18) gives a second order linear elliptic 
system in the unknowns _u. Clearly, E is the elasticity tensor and is 
assumed here to have constant coefficients (and nice "ellipticity" 
properties). Its inverse (compliance tensor) will be denoted by C. Hence: 

(2.19) x = E:e(v) < > e(v) = C:x . 

We assume again the simple boundary conditions: 

(2.20) u = 0 on 3D. 

This, of course, is strongly unrealistic: usually one has ii = ii given 

on T_. and o-n = t given on T„ . However, the proper way of 

dealing with realistic boundary conditions coincides with the one used in 
Example 1: we chose then to give more details (with simpler notations) while 

discussing the linear elliptic problem. 

One can notice here that the splitting of the problem in more than one 
unknown is natural with solid physical reasons. This is probably why the 
first mixed formulations were used in elasticity theory. We shall present 
here only one mixed formulation, which is similar to the formulation (2.9) - 
(2.10) for a single elliptic equation. We set: 

2 

H(div;D) = {t | T € (L 2 (D)) d ; T = T V i,j; div r € (L 2 (D)) d > 


and we consider the problem: 
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Find u_ e (L^(D))^ and £ £ H(dLv;D) such that 

(2.21) / (C :a ) :t dx + J u*d£V x d x = 0 ¥ £ € H(div;D), 

D - - — d — — 

(2.22) / v* div o dx = J f* v dx ¥ v € (L 2 (D)) d . 

D ~ D " 

One can see that (2.21) - (2.22) practically coincide with the varia- 
tional formulation of the Hellinger-Reissner principle. The use of this 

principle in the framework of finite elements can be traced back to the pio- 
neering work of Herrmann [33], [34], and Hellan [32]. The interest in using 

the stress field o as an independent variable is questionable in as simple 
a case as the present one, but its use is clear in more general and more 
complicated problems involving nonlinearities, plasticity and so on. 

We shall now state an abstract existence theorem that is a simplified 
version of a more general result proved in [11]. 

Theorem Is Let 5 and ¥ be real Hilbert spaces, a( ja 

bilinear form on 5x~ and b(£,i|0 a bilinear form on ~xV . Set : 

K = u| 5 € 5, b(5,<|>) = 0 ¥*€¥}, 


and assume that 


(2.23) 


3 a > 0 s.t. a (5,5) >. a»C»^ ¥ £ € K, 


b(g ,*|») 

B£H_ 


> ei+i T 


(2.24) 


3 8 > 0 s.t 


sup 

5€S-{0} 


¥ Ip e ¥ 
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Then for every € S' and €. V' there exists a unique solution 


(TJ) 

of the problem: 


(2.25) 

a(f,5) + b(5,*) = <*j,5> V 5 € 5, 


(2.26) 

b(I\i|0 = <£ 2 ,*> V * € ¥. 

■ 


Remark: Actually a stronger result Is proved in [11]. Namely: {problem 

(2.25), (2.26) has a unique solution for every € S' and e ¥'} if 

and only if {(2.24) holds and the bilinear form restricted to K, 

is nonsingular (in the sense that it induces an ismorphism from K on 
K')}. Clearly, if one assumes that a ( is symmetric and positive 
semidef inite, then (2.23) and (2.24) are necessary and sufficient for the 
existence and uniqueness of the solution of (2.25) - (2.26). 

Remark: It is clear that, if aUj,^) is symmetric, the solution 
(^jijT) of (2.25) - (2.26) minimizes the functional 

(2.27) J(5) = y a(5 ’ C) “ <£ 1’ 5> ’ 

on the subspace of 5 : 

(2.28) KU 2 ) = U\ b(C,t|0 - <i 2 A> V^ey}, 

and the formulation (2.25) - (2.26) corresponds to the introduction in (2.27) 
- (2.28) of the Lagrange multiplier ij>. 
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3. DISCRETIZING A MIXED FORMULATION 

Let us deal first with the abstract framework (2.25) - (2.26). Assume 
that we are given two sequences ^ S h^h>0 and { Vh>0 of subspaces of 
5 and V respectively. We set 


(3.1) 


*h ■ « h l «h £ 5 h- - 0 v *b e V' 


We have the following approximation theorem ([11]), 


Theorem 2: Assume that 


(3.2) 


a >0 s.t. a(5,5) > <* h »5l = V^eK, 


(3.3) 


0 h > 0 8 - fc - r ! up fm ^ifl^^h'* 1 * V * € V 

CCa^-10} a 


Then for every Jtj € 5' and € and for every h > 0 the discrete 

problem: 


(3.4) 


i(f h ,5) + b(5,* h > - <tj,5> V5€H, 


(3.5) 


b(5 h ,*) 


= <£ 2 ,*> V i|> € T h> 


has a unique solution. Moreover, there exists a constant ^h^ a h*^h^ ^ ^ 
such that 


(3.6) 


ii c - ? h » H + - * h i f i. V inf 

? h e5 h 


15 - 5 h l- + Inf I* - ^1,). 

+h € * ■ 
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The dependence of y^ on and g^ can be easily traced (see 

[11]). Clearly if (3.2) and (3.3) hold with constants a, IT independent 
of h, then (3.6) holds with a constant y" independent of h. More general 
versions of Theorem 2 (and also of Theorem 1) can be found for instance in 
Falk-Osborn [21] or in Bernardi-Canuto-Maday [9], 

We are now going to see the implications of Theorem 2 in the examples of 
the previous section. 

Example lh. Discretizations of the mixed formulations for linear elliptic 
operators. 

Many examples of successful discretizations of (2.9) - (2.10) are 

known. The first ones were introduced by Raviart and Thomas [44] and then 

elaborated and extended to more general cases by Nedelec [41], Other families 
of possible discretizations were introduced years later by Brezzi, Douglas, 
and Marini [15] and then elaborated and extended in several more recent papers 

(see, e.g., [13], [42], [14]). All of them share a very helpful property, the 

so-called "commuting diagram property", whose importance was first fully 

recognized in Douglas-Roberts [19]. Let us discuss it in a particular case: 
the BDM (Brezzi, Douglas, Marini) element of degree 2 for two-dimensional 
problems (D C IF?). Let be a regular sequence of decompositions of D 

into triangles. We assume for the sake of simplicity that ^Neu ~ ^ 

(2.2) and A(x) = 1. As a discretization of H(div;D) and L^(D) 

respectively we take 


(3.7) 


5 h = I € H(div;D); qj T € ( P 2 ) 2 VTf^}, 
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(3.8) 


* h = (v| V e L 2 (D); V| t e Pj V T€ T h >. 


Here and in the following P fc (S) (or simply P fc ) will denote the set of 
polynomials of degree < k on the set S. We consider now the discretized 
problem: 


find p, € E, and u, € 4\ 
•Hi h h h 


such that 


(3.9) 


/ Pj,* 9 dx + / u.div q_ d^c = / g n q_* ri dr V q € 5 
D D 3D ° 


(3.10) 


/ v div p,dx = / f v dx V v € ¥ . 
D ^ “ D h 


We define now an operator from (H*(D))^ into 


a h 


by: 



/ (q “ p 2 ds = 0 V e, edge in 7^, V p 2 e P 2 (e) 

e 

/ (q_ -M^q)dx = 0 ¥ T, triangle in 

/ (q_ - M^£)»rj>t b.j, die = 0 ¥ T, triangle in 


where kq- := * s t * ie cu ^^ c vanishing on 3T, and 

rot <(>: = ( — 3<f> /3x j , 3<f>/3x 2 ). We also define an operator from L^(D) 

into 4^ by: 


(3.12) J T (v - P v)p dx = 0 ¥ T, triangle in T h , ¥ Pl e P^T) 
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1 9 

Let us check now that div M^q_ = P^div iL f° r P. € (u(D)) . 


Actually, for all e we have: 


(3.13) 


/ v^div M^q d jc = / v^(M^q»ri)ds - J grad v^»M^qdx 
T 9T T 


= j v,q» nds - / grad v, • q_dx = j v, divcjdjc = J v. P, divqdx. 
31^“ T h T h T ntl 


It is also easy to check that the divergence operator is linear, continuous 
and surjective from (H*(D))^ onto L^(D). This can be summarized in the 
following diagram: 


(H^D)) 2 


l 2 (d) 


(3.14) 


It is easy to check that (3,14) implies, in particular, (3.2) and (3.3), but 
it is much more powerful than that. For instance, it implies: 


(3.15) lip - p ll < y lip - \pi 9 9 

" “ h (L 2 (D)) 2 “ 1 " **“ (L 2 (D)) 2 


(3.16) Hu - u, II , < Y, (HP “ Mu£“ o o + »u - P.uM , ), 

^ L 2 (D) " 2 “ (l 2 (d) ) h l 2 (d) 

1 2 

with Yp Y 2 are independent of h (whenever _p € (H (D)) ). In 
particular, with the choice (3.7) - (3,8) this yields: 


(3.17) 


"P “ PJI o o < Y, h 1 pll - 9 , 

" ^ (L 2 (D)) 2 “ 1 “ (H 3 (D)) 2 
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(3.18) Mu - ir M < Y h 2 (MuM 9 + II pll ). 

L (D) H (D) (H (D)) 2 

Note that (3.17) does not follow from the abstract error estimate (3.6). 

The commuting diagram has other nice properties. For instance, it allows 
a simple proof of error estimates in dual norms, as in Douglas-Roberts [20] or 
in Brezzi-Douglas-Marini [15]. Error estimates in L°° norms are also 
available: see for instance Scholz [45], [46] and Gastaldi-Nochetto [26], 

[27]. 

The most popular scheme for (2.9) - (2.10), that is, the "lowest order 
Raviart-Thomas", can be obtained by using, instead of (3.7), (3.8): 

{q| q € H(div;D); q| T € (P^T)) 2 V T; q.£| e € P Q (e) ¥ e} . 

{v | V e L 2 (D); V| t g p q v t g T h >. 

Accordingly, one then uses Pg(e) instead of ? 2 (e) in (3.11)1) and drops 
(3. ll)ii) and iii); similiarly one uses Pq(T) instead of Pj(T) in 
(3.12). It follows immediately that (3.13) still holds, and then (3.14) also 
holds. Clearly, only an 0(h) rate can now be achieved in both (3.17) and 

(3.18) . 

Example 2h. Discretizations of the Stokes equations. 

Life is much harder when we go from (2.9) - (2.10) to (2.14) - (2.15). 

The only positive aspect is that now the bilinear form a(£, v_) is such that 

1 2 

(2.23) actually holds in the whole (H^(D)) (our present E) so that 

(3.2) also holds regardless of the choice of the discretization. This might 
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partially excuse all the "Stokes-thinking" people who consider (3.3) as the 
condition for mixed methods. If one tries to discretize even the easy (2,9) - 
(2.10) with a scheme that does not satisfy (3.14), one will see that (3.2) can 
cause great difficulty. 

However, coming back to Stokes, it is true that the only condition to be 
satisfied by the discretization is (3.3), which now reads: 

/ q div v dx 

(3.20) a 8. > 0: sup - — j-q > 6 Iql 

h v£= h -{0> 'I'l - h L 2 (D)/ e h 

with, if possible, 8^ independent of H. A sufficient condition is the 
following so-called Forties trick [22]: W e have to find a linear operator 

from (H*(D))2 into such that: 

(3.21) Vv€(hJ(D)) 2 , 

(3.22) / q h div(v - M^dx =0 V v € (H*(D)) 2 , V q h £ y 

We consider one example. Let be a decomposition of D into rectan- 

gles, R, with sides parallel to the axes (the use of isoparametric elements is 
obviously also allowed, but more complicated to describe), and choose: 


(3.23) S h = {v| v € (Hj(D)) 2 ; v j r (Q 2 (R)) 2 V Re y , 

(3.24) V h = {q| q € L 2 (D); / q dx = 0; qj R £ Pj (R) ¥ R€ y . 


In (3.23) Q2(R) means the set of polynomials of degree £ 2 in each 
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variable. Let us see how to construct the operator at least for a 

1 2 

smooth _v. To deal with a general _v in (Hq(D)) is just technically 
more complicated, but the philosophy is the same. In each R we set: 

(3.25) M h2L “ _L at vertices (8 conditions); 

(3.26) J (M^v - v)ds =0 on each edge (8 conditions); 
e 

(3.27) / div(llv - v)x dsc = 0 i = 1,2 (2 conditions). 

R h- - i 

We have a total of 18 conditions (note that the dimension of Q 2 is 9). It 
is easy to check that they are independent. Let us check (3.22); that is, let 
us check that 


(3.28) / div(M^v - v)p^dx =0 V p^ € P^(R). 

R 


Clearly, (3.27) implies that (3.28) holds for p^ = and p^ = X 2 * We need 
only to check p^ = 1: 


(3.29) / div(M^v - v)dx = J (M^v - v) • ri ds = 0, 

R 3 R 


due to (3.26). We can now apply (3.6) and obtain: 


(3.30) 


II u - u^ll ^ 


Ip - 


Ph 1 ’ 2 
l 2 (d)/. 


< c h‘ 


• ( n_u. 11 3 


ipi 2 ). 


There are many other known choices available for obtaining a discretization of 
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(2.14) - (2.15) that satisfies (3.20). An almost complete list of them can be 
found in Brezzi-Fortin* together with the references. In particular, Scott 
and Vogelius [47] proved that, under minor restrictions on the decomposition 
of D into triangles, one can always use a continuous velocity field of local 
degree k and a discontinuous pressure field of local degree k-1, provided 
k >_ 4. For the low degrees, a special headache is provided by the use of 
bilinear velocities and constant pressures: its convergence has been proved 

in a variety of cases (see Johnson-Fitkaranta [38], Stenberg [48], Pitkaranta- 
Stenberg [43]) but not yet in the general case. In any event, a filtering of 
the pressure field is always required to eliminate the checkerboard modes. 
General strategies for constructing discretizations that fulfill (3.20) are 
given in Boland-Nicolaides [10] and Brezzi-Pitkaranta [16]. Modifications of 
the discrete equations that allow violation of (3.3) have received 
considerable attention in recent times. Roughly speaking, in the special case 
of equations (2.14), (2.15), they consist in substituting (2.15) with the 

"perturbed" equation 

(3.31) j q div u^ + h ^ J grad p«grad q dx =0 V q € H^(D) 

D D 

(Brezzi-Pitkaranta [16]) provided that one uses a continuous pressure field in 
the finite element discretization. A modification of (3.31) of the form: 

(3.32) J q div u dx = £ a h^ / (Aii - grad p + f )» grad q d_x 

D T€T U T 

n 


*F. Brezzi and M. Fortin, Book in preparation. 
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(a = small parameter to be adjusted; h T - diameter of T) has been introduced 
by Hughes-Balestra-Franca [35]* Note that (3.31) is simpler but has a 
consistency error of order h^ that is not present in (3.32). For more 
details and additional results in this direction, see Brezzi-Douglas [12]. 

For the use of more general boundary conditions, the basic reference is 
Verfurth [49]; see also the references contained therein. 

Additional references for the Stokes and Navier-Stokes equations can be 
found in Glowinski-Pironneau [30], Glowinski [31], Girault-Raviart [28] and 

Brezzi-Fortin^ • 

Example 3h. Discretizations of linear elasticity problems 

It is difficult, in general, to find convenient finite element discreti- 
zations for equations (2.21) - (2.22). We shall briefly indicate here four 
possible ways to proceed. The first possibility is to try to construct spaces 
that verify the commuting diagram property (as in (3.14)). This has been 
possible, up to now, only by means of composite elements : that is, each 

element is split into sub-elements and one uses trial functions that are 
polynomials in each sub-element (plus suitable continuity requirements form 
one sub-element to the other). Examples of this approach can be found in 
Johnson-Mercier [37] or in Arnold-Douglas-Gupta [4]. A second possibility is 
to give up the symmetry condition that appears in the definition of 
H(d_iy;D) and to enforce it a posteriori by means of a Lagrange multiplier. 
After discretization we deal then with stress fields having only a weak 


2 


F. 


Brezzi and M. 


Fortin, Book in preparation. 
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symmetry. This idea was first used by Fraeijs de Veubeke [23] and then modi- 
fied and analyzed by Amara-Thomas [1] and Arnold-Brezzi-Douglas [3]. A third 
possibility is to change the "auxiliary function" and use a different , 
nonsymmetric, tensor field instead of a. This will in general produce some 
trouble at I\ T (if T.. ^ 0) that can be treated with the introduction 

of an additional Lagrange multiplier on r^ gu . We re f er to Arnold-Falk [5] 
for more details on this approach. Finally, a fourth possibility is the addi- 
tion of a stabilizing term in the style of (3.31) or (3.32). For this we 
refer to Franca-Loula-Hughes-Miranda [25]. 

It should be pointed out that additional difficulties arise when dealing 
with nearly incompressible materials. In these cases (2.23) ceases to hold 
(in the limit) in the whole space but still holds for divergence-free tensor 
fields. This implies that (3.2) must also be checked if the discretization is 

such that K, <£ K. Additional references for the above (and many other) 
n 

O 

applications can be found in Brezzi-Fortin . 

4. NUMERICAL METHODS FOR SOLVING THE DISCRETIZED PROBLEM 

The major difficulty that arises in solving a linear problem such as 
(3.4) (3.5) is that the associated matrix 
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is indefinite. There are many ways of overcoming this difficulty, mostly 

using some particular feature of the problem under consideration in order to 

rewrite it in a different form. Here we shall briefly sketch two of them: 

one which is mostly used in Examples 1 and 3 and the other used in Example 2. 

The first technique, which is very old (see Fraeijs de Veubeke [24]) 

starts from the following simple observation. If the space 5, is made of 

h 

functions that are completely discontinuous from one element to the other, 

then the most natural choice of basis functions for 5^ will produce a 

matrix A in (4.1) which is block-diagonal. Then the inverse matrix A”* 

can be easily computed explicitly. Solving (3.4) element by element for 

and substituting into (3.5) (static condensation) leaves us with the final 

matrix B^A“^B and the only unkown Ij^. Note that (3.2) and (3.3) will 

imply that B T A“*B is symmetric and positive definite (if a^j,^) is 

symmetric). Now if instead the functions in 5, have some continuity 

n 

properties from one element to the other (as in Example lh we had p, *n 

— h — 

continuous at the interfaces), this cannot be done. However, one can choose 
to work in a larger space, say made of discontinuous functions, and 

then require the continuity by means of a Lagrange multiplier. Let us see the 
procedure in the particular case of Example lh. We set: 

(4.2) 5 h = (q| q £ (L 2 (D)) 2 ; q_j T £ (P 2 ) 2 V T € y , 

(4.3) = { V \ q| e € P 2 (e) V e, internal edge in y = 0 on 3D}, 

c(q_,y) - l / (jq*n)y ds. 

T€T, 3T 

n 


(4.4) 
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' Clearly, If q € a, then 

! ” h 

I 

(4.5) £ € H h C=> c(£,y) =0 ¥ y € A h . 

It is not difficult to check that the new problem 

find p e 5, , u, € v X, € A, such that 
—n h h h h h 

(4-6) / Pu*£ d * + l I «>, div 3. d £ = / go 3 .*£ d s + cCq»^ h ) v q e s. 

D TT 3D 

(4.7) I / v div P, dx = J f v dx ¥ v € 4*, 

T T D h 

(4.8) ctf^.y) =0 ¥ y € A h 

/V Arf 

has a unique solution, and that = p^, u^ = u^* Now both the unknowns 

p^ and u^ are a priori discontinuous and they can be eliminated, at the 

element level, by static condensation. The final matrix, in the unknown X^, 

will be symmetric and positive definite. It is clear that X, itself 

h 

should be an approximation of u at the interfaces, and it has been used as 

such by engineers. However, it was only rather recently that it was proved 

mathematically that X, converges to u, and in general with a better order 

h 

of convergence than u^ itself (see Arnold-Brezzi [2]). For instance, in the 
present case, once X^ is known one can construct, element by element, an 

* , V 

U h € T ) Suc h that 
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X h )p 2 ds = 0 V e edge of T, Vp 2 £ P 2 (e) 

! (u* - u h )dx = 0 
T 

and show that 

(4.10) II u - uf II , < 0(h 4 ) 

h l/(D) “ 

instead of (3.18) (for the proof of (4.10) see [15]). A similar result (but 
only with O(h^) in (4.10)) can also be achieved with the lowest order 
Raviart-Thomas element described at the end of Example lh. However, the best 
way to compute the solution for this last element is to solve with the so- 
called P^-nonconf orming method and then use the postprocessing of Marini 
[39]. For additional results on the convergence of X^ to u see Brezzi- 
Douglas-Marini [15], Brezzi-Douglas-Fortin-Marini [14] and Gastaldi-Nochetto 
[27] . 

This same idea (disconnect and use a Lagrange multiplier to force 

back the continuity) can be used for elasticity problems and in many other 

cases. However, it has not been possible, so far, to use it, for instance, 

for the Stokes equations (and more generally when continuity at the vertices 

is required in E, ). Then one can use the following other trick, that was 

h 

first analyzed by Bercovier [8]. If the space is made of discontinuous 

functions (as it was the case in our Example 2h) , then one can perturb equa- 
tion (3.5) into 



b(F h ,*) = + <£ 2 ,iJ>> v Ijj e 4' h . 


(4.11) 
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The corresponding matrix 

(for (3.4), 

(4.11)) becomes, roughly. 


/ A 

B \ 

(4.12) 

( T 

) * 


\ B 

-el 


Now the discontinuity in V allows us to eliminate \JT, at the element 

n h 

— 1 T 

level. We obtain in that way a matrix A+e BB . If (3.2) and (3.3) are 
satisfied, this new matrix will be symmetric and positive definite (always 
if a(5j,£ 2 ) is symmetric). Moreover, calling (5^,i|;^) the solution of 

(3.4) and (4.11), one has 


(4.13) I5 h - 5^1- + l* h “ - 0(e). 

The method can also be applied when ^ consists of continuous functions, 
provided that some lumping procedure is used to compute the Inner product in 
(4.11). However, in such cases, one gets for e BB a bandwidth that is 
generally larger than that of A, and this is often a considerable drawback. 

A different attempt to reduce (3.4) (3.5) to a single equation in the 

case of the Stokes problem can be found in Bramble^# 


*J. H. 


Bramble, to appear. 
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